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Abstract 

It is well known that the linear-noise approximation (LNA) exactly agrees with the chemical 
master equation, up to second-order moments, for chemical systems composed of zero and hrst- 
order reactions. Here we show that this is also a property of the LNA for a subset of chemical 
systems with second-order reactions. This agreement is independent of the number of interacting 
molecules. 
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I. INTRODUCTION 


The Markovian description of chemical systems, as given by the chemical master eqnation 
(CME) dj, is rarely exactly solvable. This has led to the development of a nnmber of 
approximation techniqnes snch as moment-closnre approximations PE] and the linear- 
noise approximation (LNA) [H |5]. The accnracy of these techniqnes is freqnently nnknown 
and hence there is a snbstantial ongoing effort to clearly nnderstand the limitations of the 
techniqnes and the magnitnde of the error inherent in their predictions (see for example 

[HHin]). 

Here we exclnsively focus on the LNA. The LNA originates from the system-size expan¬ 
sion first introduced by van Kampen [5]. The system-size expansion effectively constitutes 
an inhnite series expansion of the moments of the probability distribution solution of the 
CME in powers of the inverse volume of the compartment in which the chemical system is 
conhned. The LNA is the leading-order term in this expansion which implies that in this 
approximation, the mean concentrations are the same as given by the conventional rate equa¬ 
tions while the variance of concentration fluctuations is proportional to the inverse volume. 
Hence the LNA is conventionally regarded as a large volume (macroscopic) approximation 
of the moments of the CME; equivalently the LNA an be viewed as an accurate approxi¬ 
mation in the limit of large molecule numbers, since the macroscopic limit is the limit of 
large volumes at constant concentration |3]. Remarkably, however, the LNA is exact (up 
to second-order moments) for systems composed of zero and hrst-order reactions. This is 
since in this case the propensities (the transition rates in the CME) are already linear in the 
molecule numbers and hence the linearisation procedure inherent in the application of the 
LNA bears no effect on the equations for the time-evolution of the covariance matrix (see for 
example [H]). It has also been shown that the differences between the predictions of the rate 
equations / LNA and of the CME are proportional to the elements of the Hessian matrix 
of the rate equations |H1 El [I2| and inversely proportional to the volume. Since the Hessian 
matrix is non-zero whenever chemical systems have at least one second-order reaction, it is 
generally thought that the LNA’s predictions increase with the non-linear dependence in the 
molecule numbers in the propensities of the CME and with decreasing the volume, a claim 
supported by several case studies [HEMS]. 

Summarising, the current general picture is that (i) the LNA is exact (up to second-order 
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moments) for systems composed of at most first-order reactions since the propensities are 
constant or directly proportional to the molecule numbers, (ii) the LNA is inexact for sys¬ 
tems with at least one second-order reaction since some of the propensities are proportional 
to the product of molecule numbers. 

In this article, we show that this standard picture needs revision. In particular, we 
prove that there exists a special class of chemical systems containing at least one second- 
order reaction for which the error in LNA’s prediction of the mean concentrations and of 
the variance of fluctuations about the means is zero for all volumes and parameters of the 
system. This implies that for these systems, the LNA is exact for all molecule numbers and 
not just accurate in the limit of large molecule numbers, as commonly thought. 

The article is organised as follows. In Section II, we show by means of examples, whose 
CME can be solved exactly, that the rate equations and LNA for systems with at least a 
second-order reaction, can in some cases lead to exact expressions for the mean concentra¬ 
tions and the variance of the fluctuations about them. In Section III, we identify a broad 

class of chemical systems with this property. In Section IV, we generalise further this special 
class of systems and provide several examples of common chemical systems of this type. We 
hnish by a discussion in Section V. 

II. WHEN ARE THE RATE EQUATIONS AND THE LNA EXACT? 

A. An illustrative example 

Consider the following two different reaction systems: 

X, + X,^Xs, ( 1 ) 

ki 

Xi + X2^Xs, 0^X,, ( 2 ) 

fci fca 

where the rate constants are those which would appear in a rate equation formulation of the 
systems. Reaction ([^ is a closed heterodimerisation reaction whereby molecules of species 
Xi and X 2 reversibly combine to form a heterodimer A 3 . This system has two implicit 
conservation laws: U 2 -Ui = a and ni + n^ = f3 where Uj is the number of molecules of species 
Xi and a,P are time-independent constants. The hrst conservation law is due to the fact 
that whenever a molecule of Xi is produced (or consumed), a molecule of X 2 is also produced 
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(or consumed); the second conservation law stems from the fact that whenever a molecule 
of Xi is produced (or consumed), a molecule of X3 is consumed (or produced). Reaction 
(|^ is an open version of the previous reaction since molecules of Xi are also produced and 
destroyed. There is one implicit conservation law in this system, namely ^2 + us = 7 where 
7 is a time-independent constant. 

In steady-state conditions, both systems satisfy detailed balance and hence using standard 
methods Id the exact solution of the CME’s of both systems can be straightforwardly 
obtained. These are given by: 


F(ni) = 




fc2n (-^) 

F(ni,n2) = e ^ 


ni!(ni + a )\{(3 - ni)\M{-( 3 ,1 + a, -^) ’ 


1 + 


rii! 


ko 

( kjki \irio 

t k 2 ko ) 


hko) ’n2!(7-’^2)!’ 


(3) 


(4) 


for reaction systems ([^ and (|^ respectively. The volume of the compartment in which 
the reaction occurs is given by ff. Note that there is no explicit n2-,n^ dependence in Eq. 
(|^ since ni is related to n2,n3 via the implicit conservation law. Similarly there is no 713 
dependence in Eq. (|^ since n2 is related to 77,3 via the implicit conservation law. The 
function M{x,y,z) denotes the Kummer confluent hypergeometric function 
The rate equations for the chemical system ([^ are given by: 


A 

dt' 


!L 

dt 


02 - 


dt 


- -^00102 + 


(5) 


while for the chemical system (|^ they are given by: 


dt 

dt 


01 - ^2 - ^301 “ ^00102 + ^103) 


02 


d 

dt 


03 - -^00102 + ^103, 


( 6 ) 

(7) 


where 0 j = {ni)jfl is the concentration of species Xj. These equations when solved in steady- 
state conditions, yield the following mean number of molecules: 


{ni) = 
{ni) = 


-ako - kiQ + ’^ 40 ^ 0 ^ 1 ^ + {ako + kiQ) 


kokl 


(772) = 


2fco 

^ 1^37 


kok 2 + kikz 


(772) = ( 77 i) + a, (773) =0 -(tii), (8) 

(773) =7-(772). ( 9 ) 


for the reaction systems ([^ and (|^ respectively. 

Computing the mean number of molecules from the exact probability distribution solution 
Eq. ([^ of the CME for system ([^, one hnds that they are generally different from the 
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solution of the rate equations given by Eq. ([^. These differences are illustrated (by the 
blue and green lines) in Fig. 1(a) for the parameters ko = l,ki = 0.1, a = 0,/9 = with 
the volume 11 varied in discrete steps of 1 such that the quantity (3 is an integer (this is 
required by the conservation law discussed earlier). Note that a = 0 and /3 = H imply that 
the concentrations of Xi and X 2 are equal at all times and that sum of the concentrations 
of Xi and X 3 is unity at all times, respectively. In contrast, the mean number of molecules 
computed from the exact probability distribution solution Eq. (|^ of the CME for system 
(|^, are exactly the same as those given by the rate equation solution Eq. ([^ - this is 
illustrated by the dashed red line in Fig. 1(a). 

Similarly one can show that the Fano factors for both species (the ratio of the variance 
of fluctuations and of the mean number of molecules) computed using the exact probability 
distribution solutions differ from those obtained using the LNA for system ([^ (illustrated in 
Fig. 1(b) for the same parameters used in Fig. 1(a)) and agree exactly with those obtained 
from the LNA for system ([^. 

As mentioned in the Introduction, generally it is thought that the rate equations and LNA 
are only exact (up to second-order moments) for systems composed of at most hrst-order 
reactions and that for systems with at least a second-order reaction, they agree with the 
CME only in the limit of large volumes. From this perspective, the exact agreement found 
for system ([^, for all volumes, is perplexing since this system does possess a second-order 
reaction. Next we show that system (|^ forms part of a broad class of chemical systems for 
which the LNA is exact up to second-order moments. 


III. A SPECIAL CLASS OF SYSTEMS 


A. Specification 

Consider the following chemical system involving N species interacting via R reactions: 

fcj- 

0 —Xi, 

k- 

N k+ N 

Y,SijXi^Y^rijXi, j = 2,...,R, (10) 

i=l i=l 

where k*^~ > 0 and the following four constraints apply: 
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1. The Markov process describing the stochastic dynamics of the system is in detailed 
balance. 

2. The stoichiometric integers Sij and rij cannot simultaneously satisfy the two conditions 
rij - sij = ±1 and rij - Sij = 0 for j = 2, ...R and i = 2, 

3. The stoichiometric integers satisfy 1 < < 2,1 < ZRij ^2, j = 2, ...R. 

4. Every second-order reaction must involve Xi. 

The hrst constraint implies that the frequency of transitions from a state n to another 
state n' equals the frequency of transitions from state n' to state h, where n = {ni, n 2 ,riAr} 
and rij is the number of molecules of species Xj; note that this is not the same imposing the 
detailed balance condition on every pair of reversible reactions [1]. The second constraint 
implies that reactions j = 2,..., R cannot lead to the production or destruction of one molecule 
of Xi and simultaneously cause no change in the molecule numbers of other species. The 
latter restriction implies that reactions j = 2, ...,R cannot for example be of the type mXi ^ 
(m + l)Xi, m > 1. The third constraint implies the reactions in the system are hrst-order 
or second-order; this is since reactions involving three or more molecules are typically much 
rarer (at normal pressures and temperatures) than those involving one or two molecules. 
The fourth constraint requires that Xi must participate in every second-order reaction in 
the system, either through the binding of two molecules of Xi or else through the binding 
of one molecule of Xi and of another species. 

Note that the detailed balance condition is independent of the other three constraints 
(see later for a further discussion of this point). Hence the four constraints taken together 
imply a subset of detailed balance systems. Our claim is that for the above constrained 
chemical system, the LNA exactly agrees with the CME up to second-order moments. 

B. Statistics of molecule number fluctuations of species Xi according to the CME 

We start by showing that the steady-state fluctuations of Xi are Poissonian and uncor¬ 
related with the fluctuations of all other species, including those which participate with Xi 
in reversible unimolecular or second-order reactions. 
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It is well known that the stochastic dynamics of the chemical system (10) in well-mixed 
conditions is given by the CME [1]: 

dtP{n,t) = - l,n 2 , ■■,' 0 ,^) - klk^P^n) 

+ ki{ni + l)P{ni + 1,77.2, ••,^7v) - klniPin) 

+ n f //(h - 5,)P(h f;(n)P{n, t) 

j=2 j=2 

+ n f /r(n - Sj)Pin f-{n)P{n, t), (11) 

7=2 7=2 

where Sj is the jth colnmn vector of the stoichiometric matrix S (with elements Sij = ^ij~Sij) 
and fl is the volume of the compartment in which the chemical system is conhned. P(h, t) is 
the probability that at time t the state of the system is described by the vector n. From the 
law of mass action it follows that the (normalised) propensity functions of reaction 

j (where j = 2,..., R) are given by [1]: 


N 


//(h)=/c+ n 




N 


1 {nk - Skjy.kl^'^^ 


, fj (h) =kj n 


Hk'. 


A (uk - Tkjy . ki ^'^^ 


( 12 ) 


Now by constraint 2, the only reaction which transitions between states (tti, 77 , 2 ,..., ttjv) 
and (tti ± 1, 772 , •••, ttat) is 0 ^ Xj. Hence it follows by detailed balance (constraint 1) that 
in the CME we have the equivalence of two pairs of terms, namely: 


klkiP{ni - 1,712, = fci77,iP(77i,772, ...,7lAr), 

HfclP(77i, 77 2, ...,nAr) = fci (771 + l)P(77i + 1, 772, 


( 13 ) 


Hence the terms in Eq. (11) describing the reaction 0 ^ Xi sum to zero. The other 
terms describing reactions j = 2,..,R sum independently to zero by the detailed balance 
condition since they describe transitions other than between the states ( 771 , 772 ,..., ttat) and 


(771 ± 1,772, •••, Solving the recurrence relation given by Eq. (13), we obtain the marginal 
distribution P(ni) = X^e~^lni\ where A = klk^jk^, i.e., the fluctuations in the molecule 
numbers of species Xi are Poissonian with mean: 

VLk\ 


{ni) = 




(14) 


It also follows from Eq. (13) that 


(77i77i) = (77i)(77i) + 5i,l(77i), 7 = 1, ..., A^, 

{nlrii) = {nl){ni) + 5 i,i( 77 i)(l + 2(771)), i = 1 , N. 
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(15) 

(16) 









Thus the fluctuations of Xi are Poissonian and uncorrelated with the fluctuations of all 
other species in the system. 

We emphasize that constraint 2 is crucial to leading to Poissonian fluctuations in Xi. If for 
example we break this constraint, by considering the system of reactions 0 ^ Xi, Xi ^ 2Xi, 
then the detailed balance condition together with the CME implies that the fluctuations are 
non-Poissonian; these deviations from Poisson fluctuations are induced by the nonlinearity 
in the propensities of the reactions transitioning between states (ni,n 2 , ...,nAr) and (ui ± 
l,n 2 , However it is to be emphasised that these non-Poissonian fluctuations are also 

uncorrelated with the fluctuations of all other species in the system. Hence to summarise, 
constraint 2 leads to Poissonian and uncorrelated fluctuations in Xi while its lack leads to 
simply uncorrelated fluctuations in Xi. 


C. Statistics of molecule number fluctuations of all other species according to the 
CME 


Next we use the CME to obtain equations for the mean molecule numbers and the 
fluctuations about these means for all species besides Xi. 

By constraints 3 and 4, the reactions j = 2, ...,R in our system are first or second-order 
and involve Xi if they are second-order. This implies a simplification in the form of the 
propensities. To be specific, following are the three types of allowed chemical reactions 


and their associated propensities (for the jth reaction) as deduced using Eq. (12): (i) a 


first-order unimolecular reaction involving the decay of some species X^ is described by 
where h = 1,...,X; (ii) a second-order reaction between two molecules 
of Xi is described by ; (hi) a second-order reaction between a 

molecule of Xi and one of a different species Xh is described by (h) = . 

These three types of reactions can be generically captured by the propensity: 


..+/■ 


N 


+/- 


N 


w=2 


/-n i{ni - 1 ) ' 


(17) 


where the aX equals one if the forward reaction j is a unimolecular decay of a species X^ 
and is otherwise zero, equals one if the forward reaction j is a second-order reaction 
between two different species Xi^X^ and is otherwise zero and yj* equals one if the forward 
reaction j is a second-order reaction between the same species Xi and is otherwise zero. 





Similarly follows for the same coefficients but with minus superscript if the reverse reactions 
are of the type described. 

The time evolution of the mean concentration of species X,, (ni)/f2, is obtained by mul¬ 
tiplying the CME, Eq. 0 . by nijVt and summing over all leading to: 


(t) ^ ^ /V 

dM = ESy((/;(n)) - (/-(,!))) = 0, 

“ 3=2 


R 


= Es.jy E 

j=2 \iu=l 

R ( N 

j=2 \w^l 


a,.,3^^ + 2 ^ Pwj—F^ + Ij 


- (j^w) 


w=2 

N 

+ EA 

tt;=2 






WJ 


{nin^) ^ _ _ (ni(wi - 1)) ' 


ff2 


ff2 


(18) 


(19) 


where the angled brackets denote the statistical average. Note that in the last line we 


substituted for fj{n) from Eq. (17). There is no dependence on fcj* and since as previously 
shown, the terms in the CME describing reaction 0 ^ Xi sum to zero. The time derivative 
is set to zero since detailed balance (constraint 1) implies steady-state conditions. Using 
the fact that fluctuations in Xi are Poissonian and uncorrelated with the fluctuations of all 


other species, i.e., applying Eqs. (14)-(15), we find that Eq. (19) reduces to: 


R 


N 


N 


^ = 0 = E y + E A 


W = 1 

N 


J-2 

■E'S'yyiE 

3-2 


a 


W = 1 


WJ^ 


w=2 
N 

+ E<3 

w=2 


WJ' 




hpw + 77 01 


( 20 ) 


where we denoted the concentration of species X^, {ni)IQ, by <pi. 

Next we derive equations for the second moments of the fluctuations about the mean 
concentrations. Let Cik = kl~‘^{{nink) - {ni){nk)) be the covariance of the concentrations 
fluctuations in species Xi and X^. It can be straightforwardly shown using the latter defi¬ 
nition and Eq. 0 that the time-evolution of the covariance is given by: 


^Cik = o = y 
dt 


^(SijSkjU‘*{n) + f-{n)) Skj,. Sij 


Q 


+ + ^(K/;(h)) 


Miflin))) - (n)) - {ni){fj (h))) - (^)) - Mifj (^)))]- 


( 21 ) 


As before, this expression simplifies by the nature of the fluctuations in the number of 
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molecules of species Xi. The first sum in Eq. (21) can hence be written as: 

/ N 


R 

E 

i=2 


R 


N 


Dik{^) = X! SijSkj{fj{n) + fj (h)) = X SijSkjkH X + X l^wjMw + 7^01 

j=2 w^2 } 

R / N N \ 

(^wj(t>vo + X + Ij 4>l 

T2 


( 22 ) 




it;=2 


The main term in the second sum in Eq. (21) can be written as: 


N 


{nJlin)) - {ni){fj(n)) = kU X Kji 


{niU^) {ni){n^)' 


N 




n 


n 


X/^:. 


/ {niUin^) 

\ 

{ninj}{ni) \ {njnijni - 1 )) _ (ni(ni - l))(n^) ' 


w=2 


112 

N 


112 


112 


(23) 


W=1 


N 

I 

w=2 


- fcj'll( X Ol'^jCiw + 01 X l^wj^iw + Cii X I^wj4>w + 27 j' 0 lC'ji), (24) 


N 

E 

w=2 


where Eq. (24) follows from Eq. (23) by the application of Eqs. (14)-(16). The hrst 


sum and the last term in Eq. (24) can be easily derived. The second sum in Eq. (24) is 


less straightforward to obtain and hence we provide some additional intermediate steps, as 
follows. One hrst considers the third cumulant Ki^i which by dehnition is: 


Km = {nin^ni) - {n,„ni){ny) - {ni){nin„) - (n„)(n,ni) + 2{ni){n,j,){ni). 


(25) 


Since tc ^ 1 (as the second sum in Eq. (23) is from 2 to iV), it follows that the possible cases 


are i = l,w 1 and 1 ^ l,w ^ 1. It is easy to verify using Eqs. (15)-(16) that for each of 
these two cases, Ki^i = 0. Hence we have: 

{riiriin^) - (nin^)(n0 = {ni){{n^ni) - {ni){nj)) + {ny,){{nini) - {ni){ni)) 

= {(t)iCiyj + (t)wCii), (26) 


which leads to the second and third sums in Eq. (24) 


Hence using Eq. (22) and Eq. (24), we can deduce that the equation for the covariance 


of huctuations, Eq. (21), reduces to: 

Afc(0) ^ 


—Cik = 0 = 

dt H 


R 




(27) 


i=2 


i=2 


where A711 equals Eq. (24) and denotes the same but with minus superscripts. 


The importance of the four constraints is now clear from the derivation in this subsection 
and the previous. Constraint 1 (detailed balance) leads to fluctuations in the molecule 
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number of species Xi to be uncorrelated with the fluctuations in the molecule numbers 
of all other species. Constraint 2 leads to Poissonian fluctuations in the molecule number 
of species Xi. Constraint 3 leads to a simple quadratic form for the propensities which 
considerably simplihes the calculation. Constraint 4 leads to the equations for the second 
moments in all species to only depend on those third-order cumulants which involve Xi. By 
constraints 1, 2 and 3 we have shown that the equations for the mean molecule numbers of 


all species, Eqs. (20), are uncoupled from the second and higher-order moments; this is since 
uncorrelated fluctuations forces the condition {nin^) = {ni){nj) for ta ^ 1 while Poissonian 
fluctuations forces the condition {ni{ni - 1)) = (ni)^. By constraints 1 and 4 we have 
shown that the equations for the second-moments of the molecule numbers of all species, 


Eqs. (27), are uncoupled from the third and higher-order moments; this is since all the 
third-order cumulants involve Xi and hence must be zero since species Xi is uncorrelated 
from the rest. Thus the importance of the four constraints is that together they lead to 
the equations for the moments to naturally decouple from the higher-order moments. As 
we shall see, this property is crucial for achieving agreement of the CME with the rate 
equations and the LNA, since the rate equations depend only on the concentrations (not on 
the second-moments) while the LNA equations for the second-moments depend only on the 
concentrations and on the second-moments (not on the third-moments). 

Note that detailed balance by itself would not have led to the decoupled equations that 
we obtained, since detailed balance cannot generally guarantee Poissonian correlations in 
Xi (or in any other species for that matter - see for example Cl). and since generally the 
second-moment equations will depend on third-order cumulants which do not all involve Xi 
(see for example Appendix C of [6]). 

An important point worth mentioning regarding our derivation is that we did not use 
information about the third and higher-order moments of Xi and hence the derivation holds 
also if the fluctuations in Xi were not Poissonian but only agreed with a Poissonian up to 
second-order moments. 


D. Statistics of molecule number fluctuations using the LNA 

Next we derive equations for the same quantities using the LNA and show that they are 
one and the same as the equations we just obtained using the CME. The LNA has been 
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extensively discussed previously (see for example mm) and thus here we shall simply state 
the main results. 


As we saw earlier, constraints 1 and 2 imply that the reaction 0 ^ Xi can be treated as 
if it is separate and non-interacting with the rest of the reactions {j = 2,...,R) in chemical 


system (10). In particular the means and fluctuations of Xi can be found by considering 
the CME of reaction 0 ^ Xi while the means and fluctuations of the rest of the species are 
found by considering the CME of reactions Y,^=i SijXi ^ Y,^=i ‘TijXi where j = 2,..., R. Since 
the LNA is an approximation of the CME, it follows that we can apply the LNA to each of 
the two CMEs. 

Now it is well known that the LNA agrees with the CME up to second-order moments for 
systems of at most hrst-order reactions. Hence the application of the LNA to the CME of 
0 ^ Xi leads to the same mean number of molecules and second-moments of the fluctuations 


about this mean for species Xi as found earlier using the CME approach (see Eqs.(14)-(15)). 

Next we apply the LNA to the CME of reactions j = 2,...,R to obtain the means and 
fluctuations of the concentrations of all other species besides Xi. Given the chemical re¬ 


action system (10) and assuming that the system is deterministically monostable (see later 
for a simple proof of this property), the LNA states that the time-evolution of the mean 
concentrations is given by the conventional rate equations [1]: 




(28) 


where is the macroscopic rate vector for the forward (+) or the backward (-) reaction 

j as given by the law of mass action. Setting the time derivative to zero follows by the fact 
that constraint 1 implies steady-state conditions. Given constraints 3 and 4, the reactions j = 
2, ...,R in our system are first or second-order and involve Xi if they are second-order. This 
implies, by the law of mass action, that the macroscopic rate vector takes the following form; 
(i) a hrst-order unimolecular reaction involving the decay of some species Xh is described 
by (0) = (ph where h = l,...,iV; (ii) a second-order reaction between two molecules 
of Xi is described by = k*^~(f)'l ; (iii) a second-order reaction between a molecule of 

Xi and one of a different species Xh is described by = k'^^~<pi(j)h- These three types 

of reactions can be generically captured by the rate vector: 

( N N \ 

E + E + -it't'i . (29) 

m=l w=2 / 
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where the a’s, /3’s and 7 ’s are either 0 or 1 , depending which one of the three elementary 
interactions above describes reaction j. Note that Eq. (28) together with Eq. ( [2^ leads 
to the hnal equations determining the mean concentrations according to the LNA (and rate 
equations) and these are precisely the same as those previously obtained from the CME (see 


Eq. (@). 

As previously mentioned the LNA is only applicable if the rate equations are monostable. 
Now as we saw earlier, by constraints 1 and 2 there is only one steady-state value for the 


mean molecule number of species Xi. Furthermore Eq. (28) together with Eq. (29) are 
linear in the concentrations (f)i [i = which implies one steady-state solution for the 

concentrations of all species (monostability). 

Under the LNA, the covariance of concentration fluctuations is described by the Lyapunov 
equation [T9] : 

jQk = 0 = + ^ (30) 

at iZ 

where Jkw = Ef =2 >S'fcj^(//(0) - is the {k,w) element of the Jacobian 

matrix of the rate equations given by Eq. (28). Using Eq. (29) and the dehnition of 
above, we find 

N NR N 

WJ Ij 


E = E E E 

w=l w=l 7=2 s=2 


w=l j=2 

= t SkA 

i=2 


+/- 

ij 


(31) 


where A^j is as dehned earlier after Eq. (27). Note that here, use was made of the relation 
Cii = which is a statement of the properties of the fluctuations of Xi. 


Given Eq. (29), it is easy to verify that the hrst term in Eq. (30) is given by: 

Afc(0) = E^*A,(/;(0) + /7(0)), 

j=2 


(32) 


where Dik{4>) is as dehned in Eq. (22). 


Finally substituting Eq. (31) and Eq. (32) in Eq. (30) we obtain Eq. (27), i.e., the 
equations for the covariance of concentration huctuations according to the LNA are one and 
the same as the equations obtained earlier using the CME. 

Note that the third and higher-order moments of the concentration huctuations according 
to the LNA do not agree with those of the CME. This is since the LNA provides a Gaussian 
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approximation to the probability distribution solution of the CME [1], i.e., it leads to third 
and higher-order cumulants equal to zero, whereas the Poissonian fluctuations of the species 
interacting in second-order reactions are characterised by non-zero cumulants to all orders. 


IV. GENERALISATION AND SOME EXAMPLES 


Our proof in the previous section has been for the system of reactions (10) in which 
species Xi is special, in the sense that it is being produced and destroyed by a particular 
reaction, and participates in all second-order reactions. Relaxing this speciality of Xi leads 
to a broader class of systems for which the LNA is exact up to second-order moments. 

In particular consider the following chemical system of N species interacting via R reac¬ 
tions: 


kt kt 

fcl fcj 

N k+ N 

y ii.w. y ryA'i, 

i=l i=l 



j = S+ 1,...,R, 


(33) 


where > 0 , and R,..., *5 are positive integers taking a value between 1 and 

N. We also apply the following constraints: 


1. The Markov process describing the stochastic dynamics of the system is in detailed 
balance. 


2. The stoichiometric integers Sij and cannot simultaneously satisfy the two conditions 

“ ^ikJ = ±1 = 0 for j = S' + 1, ...R and fc = 1,..., S. 

3. The stoichiometric integers satisfy 1 < < 2,1 < Y.i'f'ij ^2, j = 2, ...R. 

4. Every second-order reaction involves at least one of the following species Xjj, ..., 


This system is a generalisation of the system (10) since the properties previously par¬ 
ticular to species Xi are now common to S different species, ..., Xj^. Constraint 

2 implies that the reactions 0 ^ Xj ^,0 ^ Xj 2,....,0 X*^ are the only reactions in the 

system which lead to the increase or decrease of one molecule of the species Xjj, Xjj, ..., 
Xjg while causing no change to the number of molecules of all other species. 
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The proof that the LNA is exact up to second-order moments for system (33) follows 
along the lines of the proof provided in the previous section. Constraints 1 and 2 lead 
to steady-state fluctuations in the concentrations of species which are 

Poissonian and uncorrelated with the fluctuations in the concentrations of all other species; 
these properties together with the special form of the propensities due to constraints 3 and 
4 lead to the exactness of the LNA in detailed balance conditions. 

A perusal of the proof provided in the previous section shows that the crucial ingredient 
for the exactness of the LNA up to second-order moments for any chemical system with up 
to second-order reactions is that the fluctuations in at least one of the species participating in 
each second-order reaction are Poissonian (up to second-order moments) and uncorrelated 


with the fluctuations of all other species; the constrained chemical systems (10) and (33) 
are examples of such systems but its unlikely that they are the only ones satisfying the 
aforementioned crucial ingredient. Furthermore it can also be deduced that if the chemical 
system only has second-order reactions between different species, i.e, = 0 m Eqs. (11) 
and (29), then exactness of the LNA is guaranteed if the fluctuations in at least one of the 
species participating in each second-order reaction are uncorrelated with the fluctuations of 
all other species (no Poissonian fluctuations restriction). 


A list of exemplary chemical systems of the type (33) and satisfying the above four 
constraints, is as follows: 

0 
0 
0 
0 
0 
0 


kn ko 

^Xi,Xi+X2^X3, 

ki ks 

(34) 

kn k 2 

^Xi,Xi + Xi^X2, 

ki ki 

(35) 

kn k 2 

^Xi,Xi + X2^2Xi, 

ki ki 

(36) 

kn ^ k 2 kA 

^ Xi, 0 ^ X 2 , Xi + X 2 ^ 2 X 2 , 

ki ki ki 

(37) 

kn k2 kA 

^X2,Xi+X2^X3^X2 + X4, 

fci ki ki 

(38) 

° Xi,Xi + X, " X,,1, * 1,...,X-1. 

(39) 


k2i+l 


These reactions describe heterodimerisation (34) and homodimerisation (35), autocatalytic 


reactions (36)-(37), an enzyme reaction (38) and a polymerisation reaction (39) leading to 
a polymer made of N monomers. 

In Appendix A, as a secondary check, we explicitly solve the steady-state CME of the 
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six reactions (34 39) in detailed balance conditions. It can be straightforwardly verified 
nsing these steady-state distribntions that the mean and variance of the molecnle nnmber 
flnctnations of all species exactly agree with those obtained from the rate eqnations and the 
LNA, respectively. From the exact solntions once can also verify that (i) the flnctnations 


of all species in the system are Poissonian and nncorrelated for chemical systems (35 37) 


and (39); this is since in each of these cases there are no implicit chemical conservation 
laws, (ii) only the flnctnations of one of the species involved in each second-order reaction is 


Poissonian and uncorrelated for systems (34) and (38); this is since these systems do have 
implicit chemical conservation laws. 

It has been shown in [20] that the network property called dehciency has implications for 
the form of the steady-state distribntion solntion of the CME and hence one might ponder 
a possible link between dehciency and the exactness of the LNA. In particnlar Anderson et 
al. showed that systems which are weakly reversible, have a dehciency of zero and lack any 
implicit conservation laws are in detailed balance and characterised by hnctnations in the 
molecnle nnmbers of all species which are Poissonian and nncorrelated, and hence by the 
resnlts of Section III it follows that for snch systems, the LNA is exact np to second-order 


moments. Reactions satisfying snch criteria are (35), (36) and (39). However generally 


systems of the type (33) and with the fonr constraints delineated above do not possess a 


dehciency of zero and they do have implicit conservation laws. For example reaction (37) 


has a dehciency of one while reactions (34) and (38) do have implicit conservation laws. 
One may also ponder whether the exact agreement np to second-order moments between 


the rate eqnations and the LNA and the CME for systems of type (33) in steady-state and 
detailed-balance conditions also extends for the whole time-evolntion of the system, i.e., can 
we relax constraint 1? Simnlations of the chemical reaction system ( [35| nsing the stochastic 
simnlation algorithm [21] conhrm that this is not the case (see Fig. 2), i.e., exact agreement 
is fonnd only in steady-state conditions. As expected, the error made by the LNA for hnite 
time decreases as the volnme of the system increases. 


A. Spatially extended systems 

Up till now we have specihcally interpreted each Xj as a nniqne chemical species, e.g. in 
the context of gene expression, Xi might refer to the mRNA whereas X 2 might refer to the 
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protein which is translated by the mRNA. Within this interpretation, the CME Eq. (0 
and the corresponding LNA given by Eq. (30) provide a non-spatial stochastic description 
of the dynamics of the well-mixed system of N unique chemical species interacting via the 


reactions (10). However by changing the manner in which one interprets Xi, one can also 
apply the results derived earlier to spatially extended systems. 

Lets say we want to model a second-order reaction A + B ^ C between a species A which 
diffuses and two non-diffusing species B and C which are localised in a certain part of space. 
This reaction could, for example, be occurring inside a cell and the localisation could be 
due to a membrane-bound compartment in which the membrane selectively allows through 
it only certain species, in this case species A 

For simplicity lets assume that space is divided into 3 equally sized subvolumes (or voxels) 
and that A is free to move between the three voxels while B and C are localised in one voxel. 
A scheme describing this spatially extended system is 0 Xi X 2 X 3 , Xi + X 4 X 5 
where Xi, X 2 and X 3 denote the same chemical species (species A) in voxels 1, 2 and 
3 respectively, while X 4 and X 5 denote species B and C respectively. Specihcally the 
aforementioned two sets of reactions respectively model the diffusion of species A from 
the surrounding space into voxel 1 and its subsequent diffusion into voxels 2 and 3, and the 
chemical interaction of species A with species B and C localised in voxel 1. 

The stochastic dynamics of this system (assuming well-mixing in each voxel) is still 
described by Eq. ( [IT| which is now referred to as the reaction-diffusion master equation 
(RDME), with the proviso that the volume H is now the volume of each voxel [23] • The 
approximate stochastic dynamics of this system is also still described by the Lyapunov Eq. 


(30) (sometimes referred to as the multi-compartment LNA [23]). Additionally this spatially- 
extended reaction system satishes the four constraints mentioned in Section III.A and hence 
it follows by the results of Section III that the multi-compartment LNA and the RDME 
agree to second-order moments in the molecule number fluctuations in each of the boxes 
composing space. It is straightforward to show that this equivalence of the two modelling 
approaches also holds if: (i) space is divided into any number of voxels; (ii) one lets species 
B and C diffuse freely throughout space; (hi) the diffusive entry of any species, from the 
surroundings into the space under consideration, occurs at any (or all) of the voxels. 


The reaction here considered is a spatially extended version of reaction (34); similarly one 


can write spatially-extended versions of reactions (35)-(39) and in each case one hnds the 
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equivalence of the multi-compartment LNA and the RDME up to second-order moments. 
By similar arguments to the one above, one can deduce that a spatially extended system 
involving N chemical species of which S of them diffuse from the surroundings into the 
space of interest, and diffuse within this space to participate in a number of reactions is a 


special case of system (33). If the chemical reactions within this space also satisfy the four 
constraints delineated previously, then the agreement of the multi-compartment LNA and 
the RDME up to second-order moments is guaranteed. 

V. DISCUSSION & CONCLUSION 


In this article, we have shown that for a class of chemical systems, the LNA gives results 
for the mean concentrations and second moments of fluctuations about the means which 
exactly match those given by the CME for all volumes, even though the chemical system 
is composed of at least one second-order reaction. In particular this also implies that the 
mean concentrations of the CME are in exact agreement with those of the conventional rate 
equations, on which the LNA is based. This is in contrast to the current prevalent thinking 
that posits the LNA is only exact for a chemical system with zero and hrst-order reactions. 
Our results are also in contrast to the fact that previous studies have found strong deviations 
between the rate equations and the CME [8l |25] and between the LNA and the CME for 
various biochemical systems involving second-order reactions dnii]. For example in [25j it 
has been found that the differences between the substrate concentrations predicted by the 
rate equations and by the CME increase as one goes further downstream in certain large 
enzyme reaction systems and in HI it was shown how the dependence of the coefficient of 
variation of protein noise versus the stress level in a gene regulatory circuit according to the 
LNA is roughly parabolic while the same obtained from the CME is a monotonic increasing 
function. What is clear from the results of this paper and the aforementioned results in 
the literature, is that the wiring of the chemical reaction network plays a major role in 
determining the deviations between the LNA and the CME, when there are second-order 
reactions. In other words, the differences between the predictions of the rate equations / 
LNA and of the CME are proportional to the elements of the Hessian matrix of the rate 
equations and also to a “wiring factor”; hence the deviations from the LNA are zero either if 
the Hessian is zero, i.e, if the chemical system is composed of zero and first-order reactions 
(the well known case) or if the system has up to second-order reactions but the wiring 
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factor is zero (the case elucidated in this paper). The four constraints on our system impose 
fairly stringent requirements on the wiring of the network and hence if they turn out to be 
necessary for the exactness of the LNA (remains to be proved) then this property is unlikely 
to be common to many realistic chemical and biochemical networks. 

Our derivation in Section III clarihes that the crucial property of the chemical systems 
for which there is an equivalence of the LNA and the CME up to second-order moments is 
that for at least one of the species participating in each second-order reaction, the steady- 
state fluctuations in the molecule numbers are Poissonian (up to second-order moments) 
and uncorrelated with the fluctuations in the molecule numbers of all species. Note that 
this property, does not exclude the possibility that one of the species participating in a 
second-order reaction (or other reaction) has non-Poissonian and correlated fluctuations. 


For example in reaction (34), there is the implicit conservation law n 2 +ns = constant which 
leads to correlation in the fluctuations of the molecule numbers of X 2 and X 3 , while in 


reaction (38), we have the conservation law ui +77,3 + 77,4 = constant which leads to correlated 
fluctuations in the molecule numbers of species Xi,X^,X 4 ^. We emphasise that for such 
cases, the exactness of the LNA extends to all species in the system, not just those exhibiting 
Poissonian and uncorrelated fluctuations. 

We have shown that the aforementioned crucial property leads to agreement between 
the CME and LNA up to second-order moments because it leads to a truncation of the 
usually inhnite hierarchy of coupled moment equations one obtains from the CME. This 
truncation is unlike that of moment-closure approximation methods Eli, in the sense that 
in the latter one imposes an ad-hoc assumption to artihcially truncate the moment equations 
whereas in our case we have shown that the truncation follows directly and automatically 
from the properties of a subset of detailed balanced chemical systems and hence is exact. 
An interesting question remains as to whether the detailed balance condition can be relaxed 
or if it is a necessary condition to guarantee the exactness of the LNA. 
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Appendix A: Detailed balance solutions of the CME of chemical systems (34 39) 
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(A6) 


The fluctuations are Poissonian for species Xi in Eq. (Al), for species Xi and X 2 in Eq. 


(A2)-(A4), for species X 2 in Eq. (A5) and for species Xi to X^r in system (A6). Note that 


the non-Poissonian marginal distributions of some species (X 2 in Eq. (Al); Xi and X 4 in 


Eq. (A5)) only originate when there is an implicit conservation law e.g. 77,2 + 773 = 77r = 


constant in Eq. (Al) and 774 + 773 + 774 = 77 ^ = constant in Eq. (A5). The exact solutions are 


obtained using the method expounded in DU. 
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Figure 1: Plots showing the dependence of the ratio of rate equation and CME mean molecule 
numbers (A) and of the ratio of the LNA and CME Fano factors (0) as a function of the volume 
for the closed hetero-dimerization reaction Xi + X 2 ^ A 3 (blue, green) and for the open reaction 
Xi + X 2 ^ A 3,0 Xi (red dashed). See text for parameter values and for the method of 
calculation. Blue and green lines and open circles denote calculations for species Xi and A 3 , 
respectively; the red dashed line denotes calculations for both Ai and A 3 . Note that lines are 
simply a guide to the eye. Since A and 0 are generally not equal to one, for the closed reaction, 
it follows that the rate equations and the LNA differ from the CME’s prediction of the first two 
moments for this reaction. In contrast A and 0 equal one for the open reaction, implying the LNA 
is exact for all volumes in this case. 
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Figure 2: Plots of the time dependence of the covariance iti ,2 in the fluctuations of species Xi and 


X 2 in reaction system (35). Solid lines show the solution of the LNA. The dots show an ensemble 
average over a large number of trajectories generated using the stochastic simulation algorithm. 
The rate constants are /cq = fei = ^2 = ^3 = 1 - The volumes are = 20 (red) and O = 40 (blue) 
in (a) and 12 = 0.3 (red) and 12 = 1 (blue) in (b). Note that the LNA agrees with the CME in 
steady-state conditions for all volumes; however the time-evolution predicted by the LNA is not 
in exact agreement with the CME and the discrepancy increases as the volume is decreased from 
12 = 40 to 0.3. 
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